WEAPONS  RESEARCH  ESTABLISHMENT  SALISBURY  (AUSTRALIA)  F/6  12/1 

an  improved  finite  element  formulation  derived  from  the  method  — ETC{U) 
JUN  77  C A FLETCHER 


30 


UNCLASSIFIED 


WRE-TM-1827(W) 


DEPARTMENT  OF  DEFENCE 


AR-000-579 


DEFENCE  SCIENCE  AND  TECJINOLOGY  ORGANISATION 


WEAPONS  RESEARCH  ESTABLISHMENT 


TECHNICAL  MEMQI 


-1827  (W)  yj 


I M ^MPROVED ^INITE  ^LEMENTJORMULATIONJ)ERIVED  FROM  j 
jjIffijilETHOD  OF  WEIGtfTEDJIESIDUALS,  j 


CM.J.yFletcher  ) x ^ O , 

ilss.®'"" 


SUMMARY 


By  replacing  the  residual  with  its  least-squares  fit  an 
improvement  has  been  made  to  the  Method  of  Weighted  Residuals. 
For  finite  element  applications  the  least-squares  fit  is  made 
on  an  element  basis.  The  improved  formulation  is  illustrated 
by  a number  of  examples  drawn  mainly  from  fluid-flow  problems. 
The  new  formulation  is  closely  related  to  the  use  of  reduced 
integration  and  detailed  comparisons  with  both  exact  and 
reduced  integration  are  presented.  The  present  formulation 
produces  superior  results  to  the  use  of  reduced  integration 
for  some  problems  and  formulations  particularly  when  used 
with  quadratic  triangular  elements.  . - 


Commonwealth  of  Australia 


June  1977 


Approved  for  Public  Release. 


CJtechnical  Memoranda  are  of  a tentative  nature,  representing  the  views  of  the 
^ ^ yuthor(s),  and  do  not  necessarily  carry  the  authority  of  this  Establishment. 


POSTAL  ADDRESS: 


The  Director,  Weapons  Research  Establishment, 
Box  21S1,  G.P.O.,  Adelaide,  South  Australia,  5001. 


w. 


A 


DOCUMENT  CONTROL  DATA  SHEET 


Security  classification  of  this  page 


1 DOCUMENT  NUMBERS 
AR 

Number;  AR-000-579 


UNCLASSIFIED 

2 SECURITY  CLASSIFICATION 


Report 
Number : 


Other 

Numbers: 


3 TITLE 


WRE-TM-1827  (W) 


a.  Complete 
Document; 

UNCLASSIFIED 

b.  Title  in 

Isolation; 

UNCLASSIFIED 

c.  Summary  in 

Isolation; 

UNCLASSIFIED 

AN  IMPROVED  FINITE  ELEMENT  FORMULATION  DERIVED  FROM  THE  METHOD 
OF  WEIGHTED  RESIDUALS. 


4 PERSONAL  AUTHOR(S): 

C.A.J.  Fletcher 


5 DOCUMENT  DATE: 


7 7.1  CORPORATE  AUTHOR(S): 

Weapons  Research  Establishment 

7.2  DOCUMENT  (WING)  SERIES 
AND  NUMBER 

Weapons  Research  and  Development  Wing 
TM-1827 

10  IMPRINT  (Publishing  establishment): 

Weapons  Research  Establishment 

1 2 RELEASE  UMITATIONS  (of  the  document); 

Approved  for  public  release. 

12.0  OVERSEAS  NO  P R.  1 A [ [ 

Security  classification  of  this  page;  I UNCLASSIFIED 


June  1977 


6 6.1 

TOTAL  NUMBER 

OFlPAGES 

> 

32 

6.2 

NUMBER  OF 

REFERENCES; 

19 

8 REFERENCE  NUMBERS 

a.  Task.  psT  76/009  ; (524) 

b.  Sponsoring 

Agency : 


9 COST  CODE: 


239  153/343 


COMPUTER  PROGRAM(S) 
(Title(s)  and  language(s)) 


’.V  . \r>  \ 


/\J 


14  I DESCRIPTORS: 

a.  EJC  Thesaurus 
Terms 


b.  Non-Thesaurus 
Terms 


Least'squares  method 
Numerical  integration 
Formulations 


Finite  element  method 
Weighted  residuals 
Reduced  Integration 


LIBRARY  LOCATION  CODES  (Cor  libraries  listed  in  the  distribution); 


SW  SR  SD  AACA  NL 


17  I SUMMARY  OR  ABSTRACT: 

(if  this  is  security  classined,  the  announcement  of  this  report  will  be  similarly  classified) 


By  replacing  the  residual  with  its  least-squares  fit  an  improvement 
has  been  made  to  the  Method  of  Weighted  Residuals.  For  finite  element 
applications  the  least-squares  fit  is  made  on  an  element  basis.  The 
improved  formulation  is  Illustrated  by  a number  of  examples  drawn  mainly 
from  fluid- flow  problems.  The  new  formulation  is  closely  related  to  the 
use  of  reduced  integration  and  detailed  comparisons  with  both  exact  and 
reduced  integration  are  presented.  The  present  formulation  produces 
superior  results  to  the  use  of  reduced  integration  for  some  problems 
and  formulations  particularly  when  used  with  quadratic  triangular  elements. 


..  ;.n  ? 


, 


Security  classification  of  this  page; 


UNCUSSIFIED 


WRE-TM-1827(W) 


1. 

2. 

3. 

4. 


5. 


1. 

2. 

3. 

4. 


5. 


6. 


1. 

2. 


3. 


4. 


5. 


6, 


7. 


TABLE  OF  CONTENTS 


Page  No. 

INTRODUCTION  1 

REDUCED  INTEGRATION  2-3 

PRESENT  FORMULATION  3-5 

EXAMPLES  5-22 

4.1  Classical  MWR  applied  to  dy/dx  - y = 0 5-7 

4.2  Finite  element  solution  of  dy/dx  - y = 0 with  one  element  7 - 8 

4.3  Two  element  solution  of  dy/dx  - y = 0 8-10 

4.4  Steady  viscous  flow  between  parallel  plates  11  - 14 

4.5  Incompressible,  inviscid  flow  14  - 22 

DISCUSSION  23  - 24 

REFERENCES  25  - 26 

LIST  OF  TABLES 

MWR  solutions  for  dy/dx  - y = 0 6 

Single  element  solutions  for  dy/dx  - y = 0 8 

Two  element  solutions  for  dy/dx  - y = 0 10 

Solution  for  viscous  flow  between  parallel  plates  13 

Origin  and  order  of  equations  (51) , (52) , (55) , (56)  and  (60) 

to  (63)  17 

Summary  of  solutions  for  inviscid,  incompressible  flow  about  a 
circular  cylinder  18 


LIST  OF  FIGURES 

Variation  of  equation  residuals  for  dy/dx  - y = 0 
Viscous  flow  between  two  parallel  plates 

Flow  field  geometry  for  two-dimensional  inviscid,  incompressible  flow 

Comparison  of  surface  velocities  - rectangular  elements  - Green’s  theorem 
applied 

Comparison  of  surface  velocities  - rectangular  elements  - Green's  theorem 
not  applied 

Comparison  of  surface  velocities  - triangular  elements  - Green's  theorem 
applied 

Comparison  of  surface  velocities  - triangular  elements  - Green's  theorem 
not  applied 


14 


1 


WRE-TM-1827(W) 


1.  INTRODUCTION 


The  main  limitations  on  the  greater  use  of  three-dimensional  finite  element 
calculations  are  excessive  CPU  time  and  excessive  space  requirements  to  achieve 
an  acceptable  accuracy.  If  it  is  possible  to  obtain  more  accurate  solutions 
from  coarse  grids  the  benefits  for  three-dimensional  calculations  are  obvious. 
For  certain  problems  reduced  numerical  integration (refs. 1,2)  produces  more 
accurate  solutions  than  exact  numerical  integration.  The  motivation  for  the 
present  work  has  come,  in  part,  from  a desire  to  explain  why  reduced  integration 
often  produces  more  accurate  solutions.  The  novel  finite  element  formulation 
presented  here  has  arisen  through  isolating  and  generalising  the  main  feature 
of  reduced  integration  that  is  responsible  for  its  success. 


2.  REDUCED  INTEGRATION 


With  the  widespread  usage  of  isoparametric  elements,  numerical  integration, 
rather  than  analytic  integration,  has  become  almost  mandatory.  Gauss  quad- 
rature formulae  have  proved  to  be  very  efficient  for  numerical  integration  over 

rectangular  elements;  an  n^^  order  Gauss  quadrature  formula  will  ensure  that 
polynomial  integrands  up  to  order  2n-l  are  integrated  exactly.  Reducing  the 
order  of  integration  introduces  an  error  into  the  evaluation  of  the  integrals 
but  it  often  produces  a more  accurate  final  solution. 

Previous  applications  of  reduced  integration  have  been  to  plate  and  shell 
problems (refs. 1,2)  to  elasto-static  problems (ref. 3) , to  plastic  flow  problems 
(ref. 4),  to  fracture  analysis(ref.5) , to  convective  transport  problems(ref .6) , 
and  to  incompressible,  inviscid  flow (ref. 7).  Reduced  integration  has  been 
compared  with  other  methods  of  directly  approximating  the  strain  field  in 
structural  applications  of  the  finite  element  method  by  Argyris  and  William(ref .8) . 
A description  of  the  general  features  of  reduced  integration  has  been  given  by 
Zienkiewicz (ref . 9) . 

An  appreciation  of  how  reduced  integration  works  may  be  obtained  by  deriving 
the  stiffness  matrix  for  a Galerkin  finite  element  formulation  in  which  the 
integration  is  performed  numerically.  Suppose  a solution  for  q is  sought  in  a 
two-dimensional  domain.  Once  an  analytic  representation,  of  the  form 


q(x,y)  = N^(x,y)  . ^,  (1) 

j 

is  substituted  into  the  governing  equation,  L(q)  = 0,  a residual  is  created. 


R(qj.x,y)  = L(q^,x,y)  . 


The  Galerkin  formulation  requires  that 


N.(x,y) 


0,  i = l,n. 


(2) 


. R(qj.x,y)  . dx  dy 


(3) 
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is  the  shape  function  corresponding  to  the  i^^  node.  Clearly  repeated 
application  of  equation  (3)  with  n different  weight  functions,  N. , will  produce 
enough  independent,  algebraic  equations  to  solve  for  the  n nodal  unknowns, 

If  the  integration  in  equation  (3)  is  performed  numerically  the  result  is 


0,  i 


1,  n. 


(4) 


th 


Wj^  is  a weight,  determined  by  the  quadrature  formula,  attached  to  the  k function 
evaluation  point.  Rj^  is  the  value  of  the  residual  at  Xj^,  yj^  and  M is  the  total 

number  of  fxmction  evaluation  points  in  the  domain.  _ 

The  residual  R.  could  depend  on  any  of  the  nodal  unknowns,  q. 

^ 3 

linear  problems. 


Thus,  for 


£ 'Sc.  ^k^  • 


j = l 


(5) 


Substitution  of  equation  (5)  into  equation  (4)  produces  the  result 
M L 


^ \ ^ ^ 'Sc.  f^k*  ^k^  • 


= 0,  i a 1,  n. 


k=l 


3*1 


(6) 


In  equations  (5)  and  (6)  L represents  the  total  number  of  nodal  parameters, 

whereas  n represents  the  number  of  nodal  parameters,  that  are  unknown. 

Equation  (5)  represents  one  linear  relationship  between  the  n unknowns,  q^ . 

Each  additional  function  evaluation  point  (Xj^,  yj^)  introduces  another  linear 

relationship  between  the  q^'s.  At  least  n evaluation  points  will  be  required 

to  ensure  that  enough  independent  linear  relationships  are  available  for  solution 
for  q^.  Directly  setting  R^  = 0 for  n values  of  (Xj^,  y^^)  would  produce  a 

solution  by  the  collocation  method,  which  is  a subgroup  of  the  method  of  weighted 
residuals.  The  Galerkin  formulation,  equation  (6),  effectively  allows  all  the 
function  evaluation  points,  M,  to  contribute  to  the  solution.  For  a solution 
to  be  possible  with  a Galerkin  formulation  it  is  necessary  that 


M > n. 


(7) 


In  practice  if  the  quadrature  formula  is  chosen  so  that  the  numerical  inte- 
gration is  carried  out  exactly  M may  exceed  n by  a factor  of  2 or  3(ref.7). 

The  excess  evaluations  of  the  residual  arise  from  the  need  for  the  solution  to 
be  constrained  by  the  local  analytic  representation  given  by  equation  (1). 

The  u.se  of  reduced  numerical  integration  implies  a smaller  number  of  evaluations, 
M,  of  the  integrand.  This  relaxes  some  of  the  constraints  on  the  solution  at 
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the  expense  of  introducing  additional  errors  by  not  performing  the  integration 
exactly.  The  empirical  evidence (refs. 1,9)  clearly  indicates  that  the  gains 
associated  with  relaxing  the  constraints  far  outweigh  the  losses  associated  with 
executing  the  numerical  integration  less  accurately. 

It  appears  from  equation  (7)  that  the  maximum  reduction  of  the  order  of  the 
quadrature  scheme  will  occur  when  M = n.  However  the  minimum  value  for  M is 
more  likely  to  be  determined  by  the  requirement  of  convergence  of  the  formulation. 
For  an  isoparametric  formulation.  Irons (ref. 10)  gives  the  minimum  order  of 
integration  as  that  which  permits  exact  integration  of  the  element  area  (in  two 
dimensions) . 

If  the  order  of  summation  in  equation  (6)  is  reversed  the  familiar  stiffness 
equation  is  obtained  i.e. 

L 

^ . q^  = 0,  i = l,n  (8) 


M 

^ "k  • "i  • "kj  '\->'k>- 

k=l 


3.  PRESENT  FORMULATION 


The  Method  of  Weighted  Residuals (ref . 11)  has  proved  to  be  a useful  framework 
for  relating  apparently  unconnected  methods  both  inside  and  outside  the  finite 
element  formulation  (ref . 12) . The  formulation  to  be  presented  here  grew  out 
of  an  attempt  to  relate  reduced  integration  to  the  framework  that  supports  the 
method  of  weighted  residuals.  Just  as  the  method  of  weighted  residuals 
includes  a broader  class  of  methods  than  finite  element  methods  so  the  current 
formulation  is  also  applicable  outside  the  finite  element  area  of  application. 

The  method  of  weighted  residuals  (MWR)  is  a suitable  technique  to  use  when 
a numerical  solution  is  sought  for 


L(v)  = 0,  (9) 

where  L is  a differential  operator.  An  approximate  solution,  u,  is  sought 
within  some  domain  D and  subject  to  boundary  conditions  on  S,  the  boundary  of  D. 
For  steady  problems,  MWR  requires  the  introduction  of  a trial  solution  of  the 
form  (in  two  dimensions) 


u(x,y) 


<^^(x,y). 


k=l 


(10) 


I 

( 

I 


rr 
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Substitution  of  equation  (10)  into  the  differential  equation,  (9),  produces 
some  residual, 

R = L(u) 


N 

■ y “k  • ‘■»k>- 

k=l 


(11) 


If  the  approximate  solution,  u,  given  by  equation  (10)  contains  the  exact 
solution,  V,  the  residual,  R,  will  be  zero  throughout  the  domain,  D.  MWR 
approximates  this  situation  by  requiring  that  the  integral  over  the  domain  of 
the  weighted  residual  is  zero.  Thus 


I 

i 

f 


R . dx  dy 


0. 


(12) 


By  repeated  evaluation  of  equation  (12),  with  different  weight  functions,  W^, 

enough  algebraic  relations  are  established  to  evaluate  the  unknown  coefficients, 
a^,  in  equation  (10).  The  present  improvement  to  MWR  consists  of  replacing  R 

by  its  least-squares  fit  over  the  domain.  Thus 


dx  dy 


0 


(13) 


is  used  instead  of  equation  (12)  to  obtain  the  algebraic  relations  between  the 
unknown  coefficients,  a^^.  R^^  is  obtained  from 


// 


dx  dy  = minimum. 


(14) 


In  one  dimension  with  <l>^  as  polynomials  in  x,  equation  (14)  can  be  differ- 
entiated to  give 


/ 


x . R . dx 


= . Rj^  . dx. 


k = 0 . . . M - 1, 


(15) 


I 


where  M is  the  order  (in  x) , of  the  residual,  R.  If  the  Galerkin  method  is 
chosen  as  the  example  of  MWR  equation  (13)  has  the  form 


R 


Is 


i 


. dx 


0. 


(16) 
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■Hius  by  setting  equation  (15)  equal  to  zero  the  replacement  of  the  residual  by 
its  least-squares  fit  may  be  interpreted,  for  this  one-dimensional  example, 
as  a generalised  Galerkin  method Cref. 12)  in  which  the  normal  Galerkin  weighting 

function,  is  replaced  by  x^  This  idea  will  be  illustrated  by  the  first 

example. 

Since  a large  number  of  efficient  finite  element  formulations  already  exist 
some  justification  for  introducing  another  would  appear  warranted.  At  the 
present  time  the  justification  is  entirely  pragmatic  : for  the  problems  con- 
sidered the  new  formulation  has  produced  more  efficient  solutions  than  a conven- 
tional Galerkin  finite  element  formulation  and  has  demonstrated  a wider  range  of 
applicability  than  reduced  integration. 


4.  EXAMPLES 


The  current  finite  element  formulation  will  be  developed  by  considering  a 
number  of  examples  of  increasing  complexity  which  will  illustrate  different 
aspects  of  the  formulation.  To  permit  an  easier  examination  of  how  MWR  works, 
and  how  it  can  be  improved,  the  first  example  will  be  solved  by  the  traditional 
MWR  formulation,  i.e.  the  approximating  functions  will  span  the  whole  domain 
rather  than  being  restricted  to  a local  region  as  in  the  finite  element  method. 

4.1  Classical  MWR  applied  to  dy/dx  - y ■ o 

A solution  is  sought,  for  the  equation  dy/dx  - y = 0 and  boundary  condition 
y = 1 at  X = 0,  in  the  domain  0 < 3^  1.  It  is  of  historical  interest  that 
this  example  was  used  by  Duncan  (ref . 13)  to  illustrate  the  traditional  Galerkin 
method. 

The  following  approximate  solution  is  introduced 

y = 1 + ai  . X + 32  . x^  , (17) 

3 


I 


I 


and  the  coefficients  aj  and  a2  are  to  be  determined.  It  may  be  noted  that 
equation  (17)  satisfies  the  boundary  condition  exactly.  Substitution  of 
the  approximate  solution,  equation  (17),  into  the  governing  equation  produces 
a residual,  R,  given  by 

R = -1  + a,  . (1  - X)  + ai  . (2x  - x^).  (18) 


Since  R is  quadratic  in  x it  is  impossible  for  ai  and  a2  to  be  chosen  so 
that  R is  identically  zero,  unless  the  choice  of  the  approximate  solution, 
equation  (17),  happens  to  contain  the  exact  solution. 

In  order  to  determine  ai  and  a2  the  integral,  over  the  domain,  of  the 
weighted  residual  is  set  equal  to  zero.  Thus 


dx  = 0,  i 


1.  2. 


If  Wj^  = x^,  i.e.  the  same  as  the  analytic  function  in 

traditional  Galerkin  method  is  obtained.  Evaluation 
to  the  following  algebraic  equations  for  ai  and  a2 , 


(19) 


equation  (1),  the 
of  equation  (19)  leads 


1 


1 


! 
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1^5  1 

6 ' ai  U ' ^2  =■  2 

1 3 1 

12  • ai  + 10  • = 3 (20) 

From  equations  (20),  aj  = 8/11  and  a2  = 10/11.  The  solution  for  y is 

Si 

shown  in  Table  1 as  y„  , . 

Gal 

TABLE  1.  MWR  SOLUTIONS  FOR  dy/dx  - y = 0 


X 

Approximate  solution 

Exact 

solution 

x 

y = e 

^Gal. 

y 

•^res,  l.s. 

^class.  l.s. 

0 

1.0000 

1.0000 

1.0000 

1.0000 

0.2 

1. 1818 

1,2057 

1.2072 

1.2214 

0.4 

1.4364 

1.4800 

1.4819 

1.4918 

0.6 

1.7636 

1.8288 

1.8241 

1.8221 

0.8 

2.1636 

2.2343 

2.2337 

2.2255 

1.0 

2,6364 

2.7143 

2.7108 

2.7183 

/ R^  dx 

0.00826 

0.00408 

0.00402 

0. 

•'o 

Substitution  of  the  solutions  for  ai  and  aj  into  equation  (18)  will 
produce  a residual  that  is  non-zero.  Intuitively  the  closer  R is  to  zero 
the  closer  y^  should  be  to  the  exact  solution.  The  conventional  way  of 

obtaining  an  improvement  in  the  solution  y is  to  allow  more  unknown 

coefficients  a^  in  equation  (17) . In  the  limit  of  having  an  infinite 

number  of  unknown  coefficients,  a^ , the  residual  can  be  made  identically 

zero  and  the  approximate  solution,  y , coincides  with  the  exact  solution. 

a 

Increasing  the  number  of  unknown  coefficients,  a^ , is  not  a very  practical 


technique  because  of  difficulties  in  solving  the  algebraic  equations  when 
the  number  of  unknowns  is  large (ref . 12) . 

A technique  for  reducing  the  size  of  R,  without  increasing  the  number 
of  unknown  coefficients,  a^ , may  be  obtained  by  approximating  the  residual, 

in  some  sense,  by  a lower  order  analytic  function.  Thus  if  equation  (18) 
is  fitted,  in  the  least-squares  sense,  by  a function  linear  in  x,  the 
result  is 


(21) 


’■4 
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Now  it  is  possible  to  ensure  that  ^ is  identically  zero  by  setting  the 


coefficients  of  x°  and  x*  equal  to  zero.  This  produces  the  result  ai  » 6/7 
and  aa  = 6/7.  The  same  result  for  ai  and  aa  would  result  from  applying  any 
of  the  methods  of  weighted  residuals  in  the  form  of  equation  (13). 


The  corresponding  solution  for  y is  plotted  in  Table  1 as  y ^ and 


it  is  apparent  that  it  is  considerably  closer  to  the  exact  solution  than 
that  produced  by  the  conventional  Galerkin  method. 

The  success  of  the  least-squares  fit  of  the  residual  is  presumably  due  to 


the  fact  that  Rj^  ^ has  the  same  global  character  as  R but  the  lower  order 


of  analytic  functions  present  has  reduced  the  constraints  on  the  unknown 


coefficients,  so  that  they  may  be  chosen  so  that  R^  ^ is  identically  zero. 


The  most  accurate  solution  might  be  expected  from  the  classical  least - 
squares  formulation,  oince  this  is  obtained  by  requiring 


1 1 
R^ 


dx  = min . 


(22) 


In  terms  of  the  method  of  weighted  residuals  this  is  equivalent  to  setting 


W.  = 3R/3a. . 

1 1 


For  the  above  problem  = (1  - x)  and  (2x  - x ). 


The 


solution  for  the  classical  least-squares  formulation  is  also  shown  in 
Table  1 and  it  is  apparent  from  both  individual  values  and  the  evaluation 


of  I R^dx  that  the  solution,  after  a least-squares  fit  of  the  residual, 
Jo 


is  very  close  to  that  obtained  from  the  classical  least-squares  formulation. 
Equation  (21)  can  be  written  in  the  form 


l.s. 


= bo  bi 


(23) 


and  the  solution  is  obtained  by  requiring  that  bo , bi  = 0.  These  conditions 
can  be  substituted  into  the  equations  that  are  used  to  calculate  bo  and  bi 
i.e. 


l.s. 


dx  = 


dx  = 0,  j =0,  1.  (24) 


1 
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Examination  of  equation  (24)  indicates  that  the  procedure  of  replac*.ig  the 
residual  by  its  least-squares  fit  is  equivalent,  for  this  problem,  to 

solving  the  original  problem  with  a modified  Galerkin  procedure  in  which  ^ 

the  normal  weight  function,  x^ , is  replaced  by  x^~^. 

4.2  Finite  element  solution  of  dy/dx  - y = 0 with  one  element 

The  equation  dy/dx  - y = 0 will  be  solved,  in  the  region  -K  x < 1 
with  the  boundary  condition  y = e“*  at  x = -1,  using  a Galerkin  finite 
element  formulation.  The  whole  domain  is  spanned  by  one  quadratic  element 

and  a solution  in  terms  of  the  equally  spaced  nodal  values,  yi  , yi  and  ya , 

is  sought  with  yi  known. 

Substituting  the  quadratic  finite  element  representation  for  y into  the 
governing  equation  produces  the  residual 


= -Jj  (x^  - 3x  + 1)  . yi  + (x^  - 2x  - 1)  . yj  - h (x^  - x - 1)  . ya . (25) 
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Application  of  a conventional  Galerkin  finite  element  formulation  i.e. 


I . R . dx  = 0,  i = 2,  3.  (26) 

•'-!  ■ 

produces  the  solution  = 1.75  yi  and  y2  = 5 yi . A least-squares  fit  of 
equation  (25)  produces  the  result 

"l.s.  = <1  >>  4>  ■ 7.  - (2X  • |)  . W M|  * i)  . (27) 

Requiring  that  R^  ^ is  identically  zero  produces  the  solution  yi  = 2.5  yi 

and  y3  = 7.yi . The  various  solutions  are  compared  with  the  exact  solution 
in  Table  2. 


TABLE  2.  SINGLE  ELEMENT  SOLUTIONS  FOR  dy/dx  - y = 0 


■ 

Approximate  solution 

Exact 

^Gal. 

y 

res.  l.s. 

X 

y = e 

-1 

0.3679 

0.3679 

0.3679 

0 

0.6438 

0.9197 

1.0000 

1 

1.8394 

2.5752 

2.7183 

If  the  integration  in  the  conventional  Galerkin  finite  element  formulation 
had  been  carried  out  numerically  a three  point  Gauss  quadrature  formula 
would  have  been  required  for  the  integration  to  be  exact.  If  reduced 
integration  had  been  used  (i.e.  a two  point  Gauss  quadrature  formula)  the 
solution  would  have  been  identical  with  that  produced  by  a least-squares 
fit  of  the  residual.  This  is  because  the  application  of  reduced  inte- 
gration to  this  problem  is  equivalent  to  replacing  the  residual  by  its 
least-squares  fit  and  either  using  the  original  weight  function  or  its 
least-squares  fit.  Since  equation  (27)  is  only  linear  in  x and  contains 
two  unknown  coefficients,  the  same  solution  will  be  obtained  whatever 
weight  function  is  used. 

4.3  Two  element  solution  of  dy/dx  - y = 0 

If  the  domain  is  not  represented  by  a single  element  the  solution  is  a 
little  more  complicated.  The  equation,  dy/dx  - y = 0,  will  be  solved  in 
the  region  0 < x < 1 subject  to  the  boundary  condition  y = 1 at  x = 0. 

The  region  is  split  into  two  elements  : element  A is  0 < x < 0.5  and 
element  B is  0.5  < x < 1.  The  solution  win  be  obtained  in  terms  of  five, 
equally  spaced,  nodal  values  yj  to  ys  with  yi  = 1.  The  shape  functions  in 
elements  A and  B are  quadratic  in  x. 

Substitution  of  the  finite  element  representation  for  y into  the  governing 
e()uation  produces  the  following  expressions  for  the  residual,  R. 
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In  element  A, 

R = (-8x^  + 22x  - 7)  . y,  + (16x^  - 40x  + 8)  . yj  (28a) 

+ (-8x^  + 18x  - 2)  . yj. 

In  element  B, 

R = (-8x^  + 30x  - 20)  . yj  + (16x^  - 56x  + 32)  . y4  (28b) 

+(-8x*  + 26x  - 13)  , y,. 


Application  of  a conventional  Galerkin  formulation  produces  the  results 

shown  in  Table  3 under  y . . 

^ ei 

To  apply  the  method  developed  in  this  paper  it  is  necessary  to  fit  the 
residuals,  in  the  least-squares  sense,  separately  in  each  element.  This 
produces  the  result: 

In  element  A, 

Rj  g = (18x  - 6.6667)  . y,  + (-32x  + 7.3333)  . yz  + (14x  - 1.6667)  . yj, 

(29a) 


In  element  B, 

R^  g = (18x  - 15.6667)  . yj  + (-32x  + 23.3333)  . y4  +(14x  - 8.3333)  . ys. 

(29b) 


By  requiring  that  Rj^  ^ in  both  elements  are  identically  zero  it  is 

possible  to  obtain  four  relationships  for  the  four  unknowns  : yz  , yj , y4 
and  ys . The  results  are  shown  in  Table  3 under  yj^^.  If  a conventional 

method  of  weighted  residual^  were  applied  to  R.  , given  by  equations  (29a) 
and  (29b),  the  results  for  y2  etc.  would  not  differ  from  y shown  in 
Table  3, 

To  produce  the  conventional  Galerkin  finite  element  solution,  shown  in 

Table  3,  by  performing  the  integration  numerically,  would  require  a three 

point  Gauss  quadrature  formula.  If  a two  point  Gauss  quadrature  formula 

(reduced  integration)  is  used  to  perform  the  numerical  integration  the 

result  is  as  shown  in  Table  3 under  y . 

r 

Examination  of  Table  3 indicates  that  the  solutions  obtained  directly 
from  requiring  Rj  ^ =0  and  from  using  reduced  integration  are  identical 

and  that  they  are  considerably  closer  to  the  exact  solution  than  is  the 
conventional  Galerkin  finite  element  solution. 

An  indication  of  the  effectiveness  of  u^ing  the_least-squares  fit  to  the 
residual  may  be  obtained  by  substituting,  y^j  and  y^^  shown  in  Table  3, 

into  the  expressions  for  the  residuals,  equations  (28)  and  (29).  The 

• I 

various  combinations  are  shown  in  Figure  1.  , R . refers  to  equations  (28). 

Jo 
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It  is  interesting  to  observe  that  generally  smaller  than 

and  is  better  distributed.  This  is  also  confirmed  by  the  values  of 

[ R*  dx  shovm  in  Table  3. 


TABLE  3.  TWO  ELEMENT  SOLUTIONS  FOR  dy/dx  - y - 0 


X 

Approximate  solution 

(nodal) , y 

Exact 

solution 

X 

y a e 

Exact 

nvimerical 

integration 

^ei 

Element 
least- 
squares 
fit  of 
residual , 

yis 

Reduced 

numerical 

integration 

^ri 

0 

1.0000 

1.0000 

1.0000 

1.0000 

0.25 

1.2707 

1.2838 

1.2838 

1.2840 

0.50 

1.6403 

1.6486 

1.6486 

1.6487 

0.75 

2.0990 

2.1165 

2.1165 

2.1170 

1.00 

2,6938 

2.7180 

2.7180 

2.7183 

0.00142 

0.00021 



0.00021 

0. 

- u - 
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4.4  Steady  viscous  flow  between  parallel  plates 

The  application  of  the  present  method  in  two  dimensions  is  not  quite  so 
straightforward.  Steady  viscous  flow  between  parallel  plates  has  been  used 
previously(ref . 12)  to  illustrate  the  Galerkin  finite  element  formulation. 

The  governing  equation  for  this  problem  can  be  reduced  to  a Poisson  equation. 
If  the  boundary  conditions  are  chosen  so  that  the  problem  has  an  exact 
solution  the  following  specification  can  be  obtained. 


d^ui  9^u  -ir/2  IT*  JT 

+ = e cos  — . 

dx^  dy^  4 2 

subject  to  the  boundary  conditions 

u = Oony  = ±1  and  x = 1 


(30) 


u = (1  - e cos  j y on  X = 0.  (31) 

The  exact  solution  is  u = cos  j y.  u is  a modified 

horizontal  velocity  difference. 

A related  problem  will  be  considered  here  because  it  leads  more  directly 
to  the  problem  considered  in  example  5.  Equation  (30)  can  be  replaced 
by  the  two  first  order  equations. 


and 


dr  ds 
dx  * dy 


(32a) 


dr  ds 
dy  dx 


0, 


(32b) 


where  r = du/dx  and  s = du/dy. 

0 < X < 1 cind  0 < y ^ 1 subject 


A solution  will  be  sought  in  the  domain 
to  the  boundary  conditions, 


-it  It 

r = . cos  j . y on  X = 0 

r = 0 on  y = 1 (33) 

s = 0 on  X = 1 and  y = 0. 

One  rectangular,  quadratic  Serendipity  element  is  introduced  to  cover  the 
domain,  as  indicated  in  Figure  2,  and  r and  s are  given  a conventional 
finite  element  representation. 


I 


1 

3 


( 


; 1 


! 
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•g-0 


FLOW  Fff  LO  MODELLED  DY  ONE  tENENOIFlTy  ELEMENT 
r ■ aw/Sx,  • • Eu/av 


Figure  2.  Viscous  flow  between  two  parallel  plates 

This  problem  has  six  nodal  unknowns  and  Galerkin  equations  based  on 
equations  (32)  are  formed  as  follows 


I ^ 3N.  — ON. 


— -jr/2  It  L . . n 

s^  - e ' • 4 • 2 ’ ^ 


i s 2,5,6 


X 3N^  _ 8 an.  _ 

E V • 'j  • E"^  ■ ° 


0,  i . 4,  7.  8. 


The  integrations  are  evaluated  numerically;  a 3 x 3 Gauss  quadrature  formula 
produces  exact  integration.  In  Table  4 are  shown  results  for  exact 
numerical  integration  and  reduced  numerical  integration  (2x2  Gauss  quad- 
rature formula) . 

The  residuals  in  equations  (34)  and  (35)  have  been  fitted  in  the  least- 
squares  sense  with  representations  of  the  form 

*^is  “*<>■*■  X3-x*y*  (36) 


Since  each  equation  is  characterised  by  four  parameters  which  depend  on  the 
six  nodal  unknowns  it  is  not  possible  to  choose  the  nodal  unknowns  in  such 
a way  that  are  identically  zero  for  both  equations  throughout  the  domain. 
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Consequently  a Galerkin  formulation  of  the  form, 

.*  .1 

/ N.  . R . dx  dy  = 0 

1 Is  ' 

'0  "^0 


(37) 


has  been  applied  to  each  equation.  The  subsequent  solution  is  shown  in 
Table  4.  It  can  be  seen  that  the  solutions  obtained  from  a least-squares 
fit  of  the  residual  and  from  the  use  of  reduced  integration  are  of  comparable 
accuracy  and  both  are  more  accurate  than  the  use  of  exact  integration. 


TABLE  4.  SOLUTION  FOR  VISCOUS  FLOW  BETWEEN  PARALLEL  PLATES 


Nodal 

parameter 

Approximate  solutions 

Exact 

solution 

Exact 

numerical 

integration 

Reduced 

numerical 

integration 

L.S.  fit 
of  the 
residual 

T2 

0.4103 

0.3698 

0.3728 

0.3265 

Ts 

0.7409 

0.6788 

0.6800 

0.7162 

0.2621 

0.2172 

0.2172 

0.2390 

S4 

1.5068 

1.4234 

1.4194 

1.2443 

Si 

0.4313 

0.4133 

0.4114 

0.3897 

St 

0.8151 

0.8100 

0.8123 

0.8798 

The  application  of  the  Galerkin  method  results  in  algebraic  equations  of 
the  form 

(K  ] (q  J = (B  1 . (38) 

where  (K  ] is  the  stiffness  matrix  and  [q  ) is  the  vector  of  nodal  unknowns. 
The  stiffness  matrices  of  the  reduced  integration  formulation  and  the  least- 
squares  fit  of  the  residual  formulation  are  identical;  however  the  matrices 
(B  ] are  slightly  different. 

A typical  term  in  the  stiffness  matrix  is  given  by 

3N. 

" J I ^i  • 

•^0  ^ 

Both  N^  and  3N./3x  are  quadratic  in  x and  y and  therefore  the  product  will  be 

integrated  exactly  by  a 3 x 3 Gauss  quadrature  formula.  If  a bilinear  least- 
squares  fit  of  3N^/3x  is  made  it  will  coincide  with  the  values  of  3N^/3x  at 

the  sampling  points  of  a 2 x 2 Gauss  quadrature  formula(ref . 14) . 

The  product  N^^  and  (dN^/dx)^^  is  bicubic  for  which  a 2 x 2 Gauss  quadrature 


formula  produces  exact  integration.  Therefore 
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,-i,-a  3N.  rV* 

Jo  Jo  h ‘ = ll  . dx  dy  . (40) 

reduced 

A contributing  term  to  [ B 1 is  the  integration  of 


/ / . cos  j y . dx  dy  . (41) 

•'0  *^0 

Clearly  a least-squares  fit  of  cos  j .y  of  the  form  given  by  equation  (36) 

will  not  necessarily  coincide  with  the  values  of  cos  j .y  at  the  sampling 

points  of  a 2 X 2 Gauss  quadrature  formula.  If  the  right  hand  side  of 
the  governing  equation  (32a)  were  a polynomial  of  second  order  or  less 
then  the  solutions  from  the  least-squares  fit  of  the  residual  and  from  the 
use  of  reduced  integration  would  be  identical. 

4.5  Incompressible,  inviscid  flow 

The  final  example  to  be  considered  is  that  of  incompressible,  inviscid 
flow  past  a two-dimensional  circular  cylinder;  the  flow-field  is  shown  in 
Figure  3. 


Figure  3.  Flow  field  geometry  for  two-dimensional  inviscid, 
incompressible  flow 

This  example  is  more  complex  than  those  considered  previously  because  the 
flow-field  is  represented  by  many  elements  and  an  isoparametric  formulation 
is  used  to  fit  the  curved  nature  of  the  body.  However  this  problem,  like 
the  previous  examples,  possesses  an  exact  solution  sc  that  direct  comparison 
of  various  formulations  is  possible.  This  problem  has  been  used  previously 
to  obtain  a systematic  comparison  of  various  elements  and  shape  functions 
for  both  exact  and  reduced  integration(refs.9,15) . 

Results  will  be  presented  here  that  compare  the  least-squares  fitting  of 
the  residual  with  exact  and  reduced  integration  for  quadratic  rectangular 
and  triangular  elements  and  two  alternative  Galerkin  formulations. 

The  governing  equations,  for  inviscid  incompressible  flow  in  two  dimensions, 
are  taken  to  be 


9u  By 
3x  * 3y 


= 0 


(42) 
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and 


9u  Jv_ 
9y  ’ 9x 


C43) 


The  Galerkin  finite  element  formulation  proceeds  by  introducing  the  following 
representation  for  u and  v 


and 


and  requiring  that 


(44a) 


(44b) 


. dx  dy 


0,  k 


1,  2 and  i = 1,  n . 


(45) 


In  equations  (44)  and  (45) , 


N. 

1 


is  the  shape  function  at  the 


.th 

1 


node. 


Uj,  are  the  nodal  values  of  the  velocity  components,  u and  v, 
and  n is  the  total  number  of  active  nodes. 

The  residuals,  Rj^,  are  obtained  by  substituting  equations  (44)  into  equations 
(42)  and  (43) , the  result  is 


and 


(46a) 


(46b) 


Substitution  of  equations  (46)  into  equations  (45)  and  rearrangement  gives 


and 


0,  i = 1.  n 


0,  i = 1,  n . 


(47) 


(48) 


A fuller  description  of  the  derivation  of  equations  (47)  and  (48)  is  given 
elsewhere (ref . 15) . 
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The  first  set  of  results  to  be  presented  have  been  obtained  with  quadratic. 
Serendipity,  rectangular  elements.  The  form  of  a^^  and  b^^  in  equations  (47) 

and  (48)  depend  on  whether  Green's  theorem  has  been  applied.  Two  cases  will 
be  considered: 

Case  1:  Green's  theorem  applied  to  the  Galerkin  equations 

^^i  . N.  . dx  dy  - 

I 3x 


a.  . 
ij 


N.  . N. 
1 J 


1 


ds 


and 


Oij 


dN. 

— ^ . N.  . dx  dy 

dy  J 


N.  . N. 
1 J 


ds . 


(49) 


(50) 


The  line  integrals  can  only  contribute  if  the  i node  lies  on  the  boundary 

and  even  then  may  produce  zero  contribution  due  to  the  boundary  conditions 

or  the  values  of  the  direction  cosines,  1 and  1 . For  an  isoparametric 

X y 

formulation  used  with  rectangular  elements  the  integrations  in  equations  (49) 
and  (50)  are  carried  out  in  the  plane  of  a regular  rectangular  element  based 
on  a i ,T}  coordinate  system.  The  result  is,  for  the  area  integrals, 


a.  . 
ij 


if 


V-'  /3N. 

3n, 

3n, 

3n,  \ 

k _ 

^ 

1 

. yk 

L \3f 

^ 1. — I 

dr} 

3t? 

3{  / 

t 

^ d^  . dJ7 


(51) 


and  b. . = 

ij 


N., 

J 


k=l 


3N.  3N,  3N.  dN, 

1 k _ 1 k 

9tj  3i7  3{ 


. X J d{  . drj. 


(52) 


X|^,  yj^  are  the  coordinates  of  the  k^^  node  in  the  element. 

Case  2:  Green's  theorem  not  applied  to  the  Galerkin  equations 

3N. 

a,  . =11  N.  . i dx  dy 


ij 


3x 


(53) 


b.  . 
ij 


3n. 

N.  . — ^ . dx  dy 
" 3y 


(54) 


For  an  isoparametric  formulation  with  rectangular  elements  equations  (53)  and 
(54)  become 


a.  . 
tj 


r / 

ir> 

L^( 

and  b.  . 
U 


N. 


3n. 

J 

3N, 

k 

3N. 

. 1 

drj 

dr} 

' H 

dN 

j 

3n, 

k 

3n. 

. _2 

3N, 

k 

dr) 

3t? 

‘ ^ 

f dt  . drj 


y di  . dn. 


(55) 


(56) 
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In  equations  (51),  (52),  (55)  and  (56)  it  is  possible  to  identify  those 
parts  of  the  equation  that  come  from  the  residual  and  those  parts  that  come 
from  the  weight  function  in  equation  (45).  Once  the  element  and  shape 
function  are  chosen  it  is  possible  to  deduce  the  order  of  the  contributions 
to  the  weight  function  and  residual.  The  results  of  this  are  shown  in 
Table  5 and  will  be  made  use  of  in  the  discussion  of  the  solutions. 


, I 


TABLE  5.  ORIGIN  AND  ORDER  OF  EQUATIONS  (51),  (52), 
(55),  (56)  and  (60)  TO  (63) 


Element 
type 


Rect- 
angular 
quad- 
ratic 
element 


Rect- 
angular! 
quad- 
ratic 
element 


Green's 

theorem 

applied 

— 

Weight  function 

Residual 

Terms  from 

Order 

of 

Terms  from 

Order 

of 

YES 

3N.  dN,  dN.  dN, 
1 ; k . i k 

dr}  3tj  3t 

bi- 

cubic 

N. 

J 

NO 

N. 

1 

bi- 

quad- 

ratic 

_ J 

3n  3n  3N  3N, 

} . .__J  . _Ji 

35  3tj  3i7  3£ 

bi- 

cubic 

YES 

3n.  3n,  3N.  3N, 

1 k 1 k 

3 L|  3 Lj  • 3 L2  3 L| 

N. 

J 

bi- 

quad- 

ratic 

NO 

N. 

1 



3n.  3n,  3N.  3N, 

3Li  3 111  SLi  3Li 

m 

The  computational  solutions,  considered  under  Example  5,  have  been 
obtained  within  the  region  ABCD  in  Figure  3.  The  nodal  points  and  the 

elements  have  been  defined  on  a polar  grid  and  an  isoparametric  formulation 
has  been  used  to  connect  this  to  a cartesian  grid.  All  results  presented 
are  for  the  variation,  with  angular  position,  of  the  tangential  velocity 
component  at  the  body  surface. 

A root  mean  square  difference,  o,  between  the  finite  element  solution 
and  the  exact  solution  at  the  body  surface  is  defined  as  follows 


o = 


- 

■ ^ 

- 

■ '“b 

- 

i=l 

- 

(57) 


where  q.p  is  the  finite  element  solution  for  the  tangential  velocity  compon- 
ent, q^  is  the  exact  solution  for  the  tangential  velocity  component  and  N^ 

is  the  number  of  nodes  between  9=0  and  90°  (B  and  C in  Figure  3).  o has 
been  found  useful  for  comparing  results  in  tabulated  form.  A summarv  of 
the  results,  for  the  various  cases  considered  under  Example  5,  is  given  in 
Table  6. 
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TABLE  6.  SUMMARY  OF  SOLUTIONS  FOR  INVISCID,  INCOMPRESSIBLE  FLOW 
ABOUT  A CIRCULAR  CYLINDER 


Element 

type 

Green's 

theorem 

applied 

Nimber 

of 

elements 

Number 

of 

unknowns 

Nodal  R. 

M.S.  differences 

Exact 

numerical 

integration 

Exact 

numerical 

integration 

Least- 
squares 
fit  of 
residual 

Quadratic 

YES 

25 

149 

0.049 

0.015 

0.040 

rectangular 

(Serendipity) 

NO 

25 

149 

0.062 

0.087 

Quadratic 

YES 

50 

199 

0.126 

0.063 

0.059 

rectangular 

NO 

50 

199 

0.252 

0.481 

0.018 

The  first  group  of  results  were  obtained  with  25  quadratic.  Serendipity, 
rectangular  elements  spanning  the  flow-field.  This  required  149  nodal 
unknowns  to  be  solved  for.  For  case  1,  in  which  Green's  theorem  is  applied, 
the  results  are  shown  in  Figure  4.  Results  have  been  obtained  for  exact 
numerical  integration,  reduced  numerical  integration  and  for  a bilinear 
least-squares  fit  of  the  residual. 
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Figure  4.  Comparison  of  surface  velocities  - rectangular  elements 
Green's  theorem  applied 
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It  is  apparent  that  the  results  for  reduced  integration  are  better  than 
those  produced  by  a least-squares  fit  of  the  residual,  and  that  both  are 
better  than  the  results  using  exact  numerical  integration.  In  obtaining 
the  solution  using  reduced  integration,  integrals  of  the  form 

r-  3N.  3N.  3n.  aN, 

(5 

Jj  J at  377  at 

are  evaluated  using  a 2 x 2 Gauss  quadrature  formula.  A 2 x 2 Gauss  quad- 
rature formula  is  capable  of  integrating  a bicubic  integrand  exactly;  the 
integrand  in  expression  (58)  is  biquintic.  Referring  to  Table  5 it  can  be 
seen  that  , which  comes  from  the  residual  is  biquadratic  and  the  term  j 

u.  u ^ L j 

which  comes  from  an  isoparametric  transformation  of  the  weight  function, 

is  bicubic. 

Sampling  at  the  second  order  Gauss  points  is  equivalent  to  replacing 
Nj  by  a bilinear  least-squares  fit  of  (ref. 14).  Therefore  sampling  the 
term  ^ in  expression  (58)  at  the  second-order  Gauss  points  may  be  inter- 


preted as  fitting  the  term  with  a biquadratic  or  bilinear  function  that 

matches  the  original  function  at  the  sampling  points. 

In  contrast  with  this  the  'least-squares  fit  of  the  residual'  results 
have  replaced  N.  in  expression  (58)  by  a least-squares  fit  but  have  left 


the  term  ^ J , that  comes  from  the  weight  function,  in  its  original  form. 

This  has  necessitated  the  use  of  a 3 x 3 Gauss  quadrature  formula.  Thus 
the  differences  in  the  reduced  integration  solution  and  the  "least-squares 
fit  of  the  residual"  solution  arise  from  the  different  treatment  of  the 
weight  function.  Why  the  reduced  integration  treatment  of  the  weight 
function  should  produce  superior  results  is  not  clear. 

For  case  2,  in  which  Green's  theorem  is  not  applied  to  the  Galerkin 
equations,  the  results  are  shown  in  Figure  5.  The  use  of  reduced  integration 
failed  to  produce  a result  and  the  least-squares  fit  of  the  residual  has 
produced  a result  that  is  inferior  to  that  produced  by  exact  numerical  inte- 
gration. 

For  this  case  it  is  necessary  to  carry  out  integrations  of  the  form 


aN.  aNj^  aN.  aNj^ 

at  ■ ■ arj  • 


Reference  to  Table  5 indicates  that  the  term 


. dt  . dr?  . 


in  expression  (59)  con- 


tributes to  the  residual.  This  term  is  bicubic  and  therefore  a least- 
squares  fit  of  this  term  will  need  to  be  biquadratic.  However  this 
requires  nine  unknown  coefficients  which  is  as  many  as  is  required  to  define 
the  overall  integration  exactly.  In  terms  of  the  required  inequality  (7)  M 
is  not  reduced.  Thus  a biquadratic  fit,  in  this  case,  violates  the  original 
requirement  of  reducing  the  number  of  constraints  on  the  residual.  The 
results  in  Figure  5 give  some  idea  of  the  error  in  replacing  the  residual 
with  its  least-squares  fit  when  no  gain  is  obtained  from  a reduction  in  the 
number  of  constraints. 
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Figure  5.  Comparison  of  surface  velocities  - rectangular  elements 
Green's  theorem  not  applied 


An  attempt  at  a reduced  integration  solution  implies  an  attempt  to  replace 
the  bicubic  residual  term  with  either  a bilinear  or  biquadratic  fit  of  the 
residual  term  whose  values  at  the  2x2  Gauss  points  match  the  original 
residual  term.  Clearly  neither  possibility  has  the  required  least-squares 
property  of  minimising  the  square  of  the  residual. 

The  second  group  of  results  were  obtained  with  50  quadratic  triangular 
elements  spanning  the  flow-field.  A solution  has  been  sought  for  199  nodal 
unknowns.  Once  an  isoparametric  transformation,  in  terms  of  the  triangular 
coordinates  Li , Lj , has  been  applied  to  equations  (49),  (50),  (S3)  and  (54) 

and  b, , are: 


the  relevant  expressions  for  a 
Case  1 : Green ' 


ij  ij 

s theorem  applied  to  the  Galerkin  equations 


a.  . 
iJ 


. !t 


j 

f / 

N.- 

) 

J J 

^ \ 

k»l 

3N. 


3N. 


dLi  3Li 


dLj  . dL2 


(60) 


and 


'ij 


> dLi  . dLj 


(61) 
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Case  2:  Green's  theorem  not  applied  to  the  Galerkin  equations 


< dN. 

1 

3n, 

k 

3Nj 

3N, 

k 

3L, 

3l, 

3 L2 

3Li 

aN. 
1 

3N, 

k 

3N. 
2 

3N, 

k 

d L] 

3 L2 

3L2 

3L, 

^ dL,  . dLj 


. > dLi  . dL2 


In  equations  (60)  to  (63)  it  is  possible  to  identify  which  parts  of  the 
equations  come  from  the  weight  function  and  which  parts  come  from  the 
residual  in  equation  (45).  This  information  is  given  in  Table  5 along  with 
the  order  of  the  various  terms  in  equations  (60)  to  (63) . 

For  case  1 above  the  results  are  shown  in  Figure  6.  Results  are 
presented  for  exact  numerical  integration,  for  reduced  numerical  integration 
and  for  a least-squares  fit  of  the  residual  of  the  following  form 

R,  = ao  + ai  • Lj  + 32  . L2.  (64) 


A seven  point  formula  was  used  to  produce  the  exact  numerical  integration 
results  and  a four  point  formula,  due  to  Cowper (ref . 16) , was  used  to  produce 
the  reduced  integration  results;  these  quadrature  formulae  are  described 
elsewhere (ref . 7) . 
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Figure  6.  Comparison  of  surface  velocities  - triangular  elements  - 
Green's  theorem  applied 
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The  "least -squares  fit  of  the  residual"  results  are  better  than  the 
reduced  integration  results  and  both  are  better  than  the  exact  integration 
results  (see  Table  6 also).  It  is  interesting  that  the  reduced  integration 
results  have  required  four  evaluations  of  the  residual  per  element  whereas 
the  "least-squares  fit  of  the  residual"  results  have  required  only  three 
parameters  per  element.  Thus,  in  terms  of  the  inequality  (7),  the  number 
of  unknowns  n is  199  for  both  cases.  The  total  number  of  function  evalua- 
tions, M,  for  reduced  integration  is  400  and  for  the  least-squares  fit  of 
the  residual  is  equivalent  to  300.  Examination  of  Table  5 indicates  that 
in  equation  (60)  and  (61)  is  biquadratic  so  that  a least-squares  fit  that 

is  linear  in  Li  and  Lj  would  appear  appropriate. 

The  results  for  case  2,  when  Green's  theorem  is  not  applied,  are  shown 
in  Figure  7.  For  this  case  the  results  using  exact  integration  are  poor 
and  the  results  using  reduced  integration  are  worse.  But  the  results 
using  a least-squares  fit  of  the  residual  are  very  good,  virtually  as  good 
as  the  best  reduced  integration  results  obtained  for  this  probem  (Figure  4) . 


ANGULAR  POSITION,  tl  (OfGREtS) 

Figure  7.  Comparison  of  surface  velocities  - triangular  elements  - 
Green's  theorem  not  applied 
The  least-squares  fit  of  the  residual  was  of  the  form 

= ao  + ai  . Li  + a2  . Lj  + a3  . Lj  . Lj . (65) 


The  term  in  equations  (62)  and  (63)  that  contributes  to  the  residual  is 
listed  in  Table  5 and  is  quadratic  in  Lj  and  Lj . 

The  relatively  poor  showing  of  reduced  integration  applied  to  triangular 
elements  may  be  attributed  to  the  fact  that  the  lower  order  numerical  inte- 
gration formula  used  has  no  least-squares  interpretation  as  has  the  corres- 
ponding Gaussian  quadrature  formula(ref . 14)  used  with  rectangular  elements. 
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5.  DISCUSSION 


Least-squares  fitting  has  beer,  used  previously  in  finite  element  formulation. 

In  the  area  of  plate  bending.  Irons  and  Tazzaque(ref.l7,18)  have  started  with 
9 and  12  degree  of  freedom  triangular  elements  and  replaced  the  second 
derivative  of  the  shape  function  by  its  least-squares  fit.  This  has  allowed 
the  use  of  a lower  order  numerical  integration  formula  and  produced  superior 
results  to  the  elements  they  started  with. 

Hinton  and  Campbell (ref . 19)  have  used  local  and  global  smoothing  in  order 
to  improve  stresses  obtained  from  numerically  integrated,  two  dimensional, 
isoparametric  elements.  Hinton  and  Campbell  note  that,  for  rectangular, 
quadratic  elements  the  evaluation  of  the  stiffness  integral 

. D . B , dA  (66) 

by  reduced  integration  (2x2  Gauss  quadrature  formula)  produces  the  same 
result  as  performing  the  exact  integration  of 

B ^ . D . B . dA  (67) 

where  ~ indicates  a local  least-squares  fit  of  that  term.  Since  each  least- 
squares  fit  is  bilinear,  the  integrand  of  expression  (67)  is  bicubic  and  can 
be  integrated  exactly  by  a 2 x 2 Gauss  quadrature  formula.  The  integrand  is 
of  the  same  order  as  that  produced  when  expression  (58)  is  replaced  by 


ff  N.  < 

■ 

3N, 
k 

0N. 
1 

0N, 

k 

11  ’1=  1 

01? 

0TJ 

(68) 


It  is  apparent  from  previous  applications,  and  from  the  results  presented 
in  this  paper,  that  reduced  integration  has  produced  superior  results  to  the 
use  of  exact  integration  for  rectangular  elements  but  has  not  been  effective 
for  triangular  elements. 

An  important  step  in  the  formulation,  presented  here  as  an  alternative 
to  reduced  integration,  has  been  to  recognise  that  if  the  residual  is 
written 


R = 


Rj  (x.y)  . qj 


(69) 


where  q^  are  the  nodal  parameters  and  Rj(x,y)  are  made  up  of  shape  function 
derivatives,  etc,  determined  by  the  governing  equations,  then  it  follows  that 

L 


R 


Is 


R. 

J 


Is 


{x,y)  . . 


(70) 


J = 1 


Equation  (70)  is  particularly  useful  since  it  is  straightforward  to  form  R. 


•hs 


on  an  element  basis.  If  the  numerical  integrations  are  performed  on  a dummy 
element,  once  and  for  all(ref.l5),  the  formation  and  evaulation  of  R.  is 


1 


also  economical. 


Is 
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It  has  been  demonstrated  that  the  use  of  a modified  method  of  weighted 
residuals 


T 


W . R,  . dA  = 0 

Is 


(71) 


has  produced  more  accurate  results  than  the  use  of 


f 


W . R . dA  = 0 


(72) 


in  a larger  number  of  situations  than  has  the  use  of  reduced  integration 
applied  to  equation  (72).  In  particular  the  current  formulation  has  extended 
the  benefits  associated  with  reduced  integration  to  triangular  elements. 

It  seems  likely  that,  where  reduced  integration  produces  very  accurate 
results,  the  increased  accuracy  over  the  use  of  equation  (71)  comes  from  the 
implicit  treatment  of  the  weight  function,  W.  This  would  appear  to  be  a 
fruitful  area  for  future  research. 
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